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1 Introduction 



Motivation Since Loup Verlet's landmark paper in 1967, the classical Verlet 
algorithm has been the main workhorse for constant energy molecular dynam- 



ics [ 53 1 . It is an attractive algorithm for the Hamiltonian ODEs that arise in this 
context because of its explicit, time-reversible, and symplectic nature. In partic- 
ular, symplecticity implies long time stability of the Verlet algorithm. The usual 
proof of this statement uses backward error analysis to show that level sets of a 



nearby Hamiltonian function interpolate Verlet trajectories [|4}[36J[37J . This prop- 
erty implies that a Verlet trajectory is confined to these level sets for the duration 
of the simulation. As a consequence a Verlet integrator nearly preserves the true 
energy and exhibits linear growth in global error. Versions of Verlet to constrained 
(RATTLE) and multiscale (RESPA) Hamiltonian systems are also available. For 
these reasons Verlet integrators are well-suited for long time simulation of con- 
stant energy molecular dynamics. 

The situation is quite different in the context of constant temperature molec- 
ular dynamics. A molecular system at constant temperature visits every energy 
isosurface with nonzero probability and its evolution is typically modeled using 
ergodic stochastic differential equations (SDE). In contrast to their deterministic 
counterpart, explicit integrators for SDEs diverge from the equilibrium behavior 
of the true solution Q. This divergence is easy to understand since explicit in- 
tegrators are only conditionally stable. Indeed for any time-stepsize one can find 
an energy above which an explicit integrator is unstable. As a result stochastic 
effects necessarily induce instabilities by driving trajectories to these energy val- 
ues. Since molecular simulations often involve unbounded potential energy (e.g., 
Lennard- Jones interaction), these energy values are attainable, and this issue calls 
for new integration strategies for constant temperature molecular dynamics. 

Constant Temperature Molecular Dynamics We briefly recall what it means 
for a molecular system to be at constant temperature. Consider n molecules with 
masses m ; for i = 1, n evolving in a d-dimensional periodic box (or torus in 
v = dn dimensions TP). Let M represent the diagonal mass matrix of the molec- 
ular system. We assume the particle interaction is given by a potential energy 
function U : T v R. The Hamiltonian H : T u x W R of this system can be 
written as: 

H(q,p) = ±p T M- 1 p + U(q), (1) 
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where q E T u and pel" represent respectively the positions and momenta of 
the molecules. In terms of this Hamiltonian, define the probability distribution fi: 

fi(dq, dp) = Z- l e- pH ^ p) dqdp Z = [ e -^ H ^ p) dqdp . (2) 

JT"xW 

Here we have introduced the parameter (3 which is inversely related to the temper- 
ature T and Boltzmann constant ks via (3 = 1/(/cbT). 

A molecular system with Hamiltonian ([T]) is at constant temperature T if its 
trajectories sample from the probability distribution (|2]). The standard way to guar- 
antee that this is indeed the case is to assume that the molecular system follows a 
continuous, stochastic dynamics of the form: 

dQ , 

— = M P 

dt ' (3) 

dP = -VU(Q)dt + drj(Q, P) , 

where drj : T u x K" — > W represents a thermostat force. Physically one can 
interpret the thermostat force as modeling interaction between the molecular sys- 
tem and a heat bath. Mathematically it is essential that the thermostat force 
be stochastic to ensure that the dynamics of ([3]) be ergodic with respect to (|2]). 
Several specific forms of the thermostat force have been proposed in the litera- 
ture [10|TTJ24 1 41 43 1 and which one is best remains open for debate. This is a 



modeling question which is beyond the scope of the present paper. Here we as- 
sume that the thermostat force is given, and we propose an integration strategy that 
does not make strong assumptions on its precise form. In the applications and the- 
ory sections of this paper we focus on drj given by Langevin dynamics fl0||4"3 | and 
in a companion paper [[8j consider other thermostats including stochastic rescaling 



dynamics pTp2| and Nose-Hoover-Langevin dynamics p4pT| . The assumption 
of continuity excludes, e.g., the Andersen thermostat because it involves discrete 
collisions at random times that in their wake leave the momentum of the molec- 



ular system discontinuous PJ. We refer the reader to [ 15 29 1 for recent progress 
quantifying the mixing properties of the Andersen thermostat for molecular sys- 
tems. 

The SDE ([3]) is characterized by degenerate noise, irregular drift, high-dimen- 
sionality {y is typically very big), and non-well-separated time-scales. In this 
context the main aim of numerical methods is to estimate long-time dynamical 
properties. This calculation is typically done by launching a single run of an 
explicit integrator and collecting statistics. However, without the patch introduced 
below this approach is prone to failure due to numerical instabilities. 



3 



To be concrete consider computing the time-correlation in momentum along 
an equilibrium path of the SDE. This computation is common in molecular dy- 
namics and the reader is referred to [|2}[T8j for expository accounts. Define the 
continuous equilibrium correlation in momentum as: 

A(r) = (P(r + t) T P(t)) , (t>0), (r>0), (4) 

where the angle brackets denote a double average with respect to realizations of 
([3]) and an initial condition distributed according to ([2]). The usual way to estimate 
A(t) over a time interval [0, T] is by a sample average computed on-the-fly using 
a single run of an integrator. The numerical equilibrium correlation is defined as 
the limit as the number of samples tends to oo: 



^ (r) = iToo [jf ^ P W)AJ P LV»J I • (5) 

The difficulty is that this limit does not generally exist if the integrator is explicit. 
This divergence is an established problem with explicit discretizations of SDEs 
that possess drifts of limited regularity 




Proposed Integration Strategy This paper proposes a new integration strat- 
egy to solve the SDE (|3]) based on combining an explicit integrator with Monte 



Carlo methods to sample from the SDE's equilibrium distribution [19 30 34]. 
The resulting 'Metropolized integrator' preserves the equilibrium distribution and 
is often provably ergodic. This feature motivates their use as sampling meth- 



ods [JT ,13, 39 , 42 , 44 1 . In addition, in [|9j we showed that a Metropolized integrator 



also approximates pathwise the SDE's solution on finite-time intervals. 

These properties ensure that a Metropolized integrator can be used to esti- 
mate dynamics along an infinitely long solution of ([3]). Indeed one can generate a 
long time trajectory of a Metropolized integrator, and along any finite-time inter- 
val update sample averages of dynamic quantities. These averages converge as a 
consequence of ergodicity of a Metropolized integrator. The averages can also be 
made arbitrarily close to the true solution's average by selecting the time-stepsize 
small enough. For example, a Metropolized integrator can be used to approximate 
to arbitrary precision the equilibrium correlation function A(r). In fact, we show 
in this paper for every T > 0, there exists a C (T) > such that for h sufficiently 
small 

sup \A h (r) - A(h[r/h\)\ < C(T)h . (6) 

re[0,T] 
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The constant C(T) increases monotonically with the length of the time-interval 
T. Hence, one cannot use this integrator in situations where T is very large like 
rare event simulation. For such problems the reader is referred to methods adapted 
to molecular systems with rare events such as milestoning (5T 52 1. 



The error estimate ([6]) provides a theoretical order of accuracy of a Metropolized 
integrator. However, for the strategy to be practical, several questions remain: 

• what does the patch involve? 

• is the patch scalable with respect to system size? 



does the patch work with RATTLE (Verlet with constraints) [ 50 1 or RESPA 
(Verlet with multiple time-step-sizes) [ |28l[48| ? 



The aim of this paper is to answer these questions. The paper is organized as 
follows: 

§|2] provides step-by- step instructions on how to patch explicit integrators based 
on Verlet, RATTLE, and RESPA, and some basic theory explaining why the 
patch works; 

§|3] conducts numerical tests on a Lennard- Jones fluid and 'dumbbell' systems; 
§|4] contains some conclusions and future improvements. 

2 Patch 

The patch involves splicing an explicit integrator with Metropolis steps. We will 
show that the resulting integrator possesses the following properties: 

(PI) preservation of the SDE's equilibrium distribution; and, 

(P2) pathwise accuracy on finite time-intervals. 

To illustrate how the patch works, we shall implement it on an explicit integrator 
based on splitting the SDE ([3]) into Hamilton's equations: 

= M l P , 

(7) 

= -VU(Q) , 
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and equations describing the effect of the thermostat: 




(8) 

= drj . 

The explicit integrator considered is defined as a composition of a step of Verlet 
for (|7]) and a step of an approximation to ([8]) (or vice versa). If Verlet is replaced 
by an implicit method (e.g., implicit Euler), then the splitting will satisfy prop- 
erty (P2). However, due to discretization error, even an implicit integrator will 



generally fail to satisfy property (PI) [45 -47 1. 



Before we continue let us introduce some notation. Let = T u x R. u denote 
the 2z/-dimensional phase space of the molecular system and Y(t) = (Q(t), P(t)) 
{} denote the true solution of ([3]) at time t > with initial condition Y(0) = x E 
n. In what follows we take for granted that this solution exists for all time. 



2.1 Explicit Integrator 

Given a time-stepsize h and time interval T, set the number of steps to be N = 
[T/h\ and introduce an evenly-spaced mesh in time tk = hk for all < k < N. 
Let ipt k+1 ,t k '■ ^ — > ^ denote an approximation to the thermostat dynamics ([8]). 
This map depends on time because of the stochastic effects in the thermostat. For 
example, for the Langevin dynamics the thermostat force is given by, 

dr] = -^M^Pdt + yj2 1 p~ l dW , 

where W is a //-dimensional Wiener process and 7 is a thermostat parameter 
[10j43|. In this case ([8]) are Ornstein-Uhlenbeck equations in momentum whose 
pathwise unique flow is almost surely: 

VW : ^ (9>e^ M ~ lfc P + *fwJ ' (9) 

where we have introduced the random vector: 



1 I k+1 e -'<M-Ht k+1 -s) dw ^ _ 
Jt k 



The map 4>t h+1 ,t k satisfies property (PI). For other SDE-based thermostats an 
approximate map that satisfies (PI) can be similarly constructed. Notice that this 
map does not alter the positions of the molecular system. 



6 



We introduce a second map 9 h : — > which approximates Q. Since the 
Hamiltonian is time-independent, this map depends only on the time-stepsize h = 
tk+i — tk. The patch we introduce below will require that this map is symmetric 
and volume-preserving [20 25 1. An explicit integrator for (|3]) is then given by: 



0t k+1 ,t k = VViA ° °h ■ (10) 



The order in ( fT0[ ) does not matter since the integrator's single step accuracy is 
0(h 2 ) either way. Higher-order accurate or implicit schemes can also be patched, 
but such integrators may require more computational effort per step. To ensure 



scalability the map 9 h in ( fTO] ) will use Verlet to separately update sets of particles 
of the molecular system. 

To this end we partition the molecular system into m sets of particles so that 
each set has Uj degrees of freedom and J2T=i u j = v - Given a time stepsize h and 
initial condition x 6 f2, a single step of ( fTO] ) is defined as: 

-Xa = i>t lt t ° Qh, m ° ■ ■ ■ ° dh,i{x) (11) 

S v " 

0,-, 

This update gives a numerical approximation to Y(h), X 1 m Y(h). The map 
0h,j is defined as a Verlet update of the position and momentum of the jth set of 
molecules fixing the other sets. More precisely given an input 

X = ((«!,••• ,Q m ),(Pir-- ,P m )) • 

where (q^Pj) G T"* x W*, the map 9 h j outputs 

8h,j(x ) = ((<Zi,-- - ,Qj-uQj,Qj+ir" ,«m)>(Pir- - ,Pj-i,Pj,P j+ i, • • • ,P m )) • 



Here: 




i- 



= ^ + hMj 1 pj i , (12) 

where JVfj is the subset of the global mass matrix M associated to the jth set of 
particles. 

The map 9hj is symmetric and volume -preserving on f2 since it is a Verlet 
update with respect to the Hamiltonian Q restricted to: 

{Qi} x • • • x x (T"0 x {q j+1 } x • • • x {qr m }x 

x ■ ■ ■ x {p._ x } x (R"i) x {p i+1 } x ■ ■ • x {pj C A . 
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Even though the composite map h is a first-order splitting of ([7]), the one-step 
error of 9hj in preserving energy is 0(h 3 ) because each of the separate Verlet 
updates is second- order accurate. 



2.2 Metropolized Integrator 

To derive an integrator that satisfies property (PI), we patch the explicit integrator 
( [TTj i with Metropolis accept or reject steps as follows. Given the time-stepsize h, 
the initial condition x G f2, and the uniform random numbers Q ~ £7(0, 1) for 
1 < j < m, one step of a Metropolized Verlet integrator determines X 1 rj Y(h) 
using: 

X! = 1p tl ,t ° 0h,m° ■ ■ ■ ° 0h,l(x) . (13) 
v v ' 

Here we have introduced the stochastic maps 9hj ■ ^ — > ^ which are Metropolized 
versions of the deterministic updates 9 h j. Given input 

x = ((qi,--- ,Q m ),(Pir-- ,P m )) ' 
the map 9hj computes a proposed move 

using a Verlet update 

P* i = Pj - f Vg^f/^!, • • • , q m ) , 

= q 3 + hMfp*! , (14) 
Pj = P*i ~ IVq^.C/^x, ■ ■ • , q^x, <? i+1 , • • • , q J , 

and accepts this proposed move with probability 

lAexp(-f3[H(x^)-H(x )]) . (15) 

If this proposal is rejected the momentum is reversed. To summarize 

9 hJ ~. xq i y X\ , 

with output X\ defined as 
x i = ((<Zi,--- ,Qj-i,Qj,Qj+i,--- ,q m ),(Pir-- ,Pj-i,Pj,Pj+i,--- ,P m )) ■ 
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where 




; , if 0<lAoq)(-/?[F(a!*)- J ff( a j)]), 

{Qj,Pj),= i / \ • ( lf) ) 

i „ „ i otherwise . 

The patch we propose consists of an accept or reject step like ( fT6] ). It requires 
evaluating the total energy at the current step and at the proposed moves. These 
statistics are usually computed alongside evaluations of the force field. When a 
proposed move is rejected, the momentum of the jth set of particles is reversed. 
While these rejections ensure the integrator is unconditionally stable and preserves 
the SDE's equilibrium distribution, they cause an 0(1) error in accuracy due to 
momentum reversals. Next we consider the effect of these rejections and momen- 
tum reversals on the approximation to the dynamics of the SDE. 



2.3 Quantitative Error Estimates 

Consider once more the Hamiltonian of the molecular system: 

H(q,p) = ^p T M~ 1 p + U(q) (17) 

where q E T u and p E W represent a configuration and momentum of the molec- 
ular system, respectively. For simplicity, we will assume the thermostat dynamics 
is given by Langevin [[T0}|43j: 

I dt . , (18) 

]dP = -VU(Q)dt - jM^Pdt + ^2-ff3~ 1 dW , 

where 7 is a thermostat parameter, (3 is the 'inverse temperature' entering the 
distribution ([2]), and W is a //-dimensional Wiener process. The theory below 
will rely on the following regularity of the potential energy. 

Assumption 2.1. The potential energy U :T"->lis smooth. 

Remark 2.2. This assumption does not permit the potential force to have singulari- 
ties. For example, it holds for Morse potential interactions, but not Lennard- Jones 
interactions. One can relax this requirement by replacing smoothness of U by 
some coercivity. 
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Let W denote expectation conditioned on the initial distribution being the 
equilibrium distribution of the SDE ( |T8| ): 



E»(g(X k ))= / E x {g{X k )) fJl (dx), X = x e ft 



n 



The following theorem states that the Metropolized integrator satisfies property 
(P2). 



Theorem 2.3. Assume \2J\ For every T > 0, there exist positive constants h c and 
C (T) such that, 

(E"{\X it/hi -Y(h[t/h\)\ 2 }y /2 < C(T)h, 

for all h < h c and t e [0, T]. 

A proof of this result relies on single-step accuracy and bounds on moments 
of the Metropolized integrator (see [9]). These bounds help to boost the single- 
step error estimate to a global error estimate, and hence, pathwise convergence 
on finite time-intervals. The Metropolized integrator initiated from equilibrium 
satisfies such bounds as a consequence of property (PI). 



Theorem 2.3 does not require that the Metropolized integrator be ergodic, but 
only that it preserves the equilibrium distribution. For this reason one can extend 
this result to Metropolis-adjusted discretizations of other SDE-based thermostats. 
Ergodicity is technically difficult to establish when the thermostat force is a non- 
linear function of the state or highly degenerate (see, e.g., [24]). However, for 



Langevin dynamics (18) with linear friction and additive noise on all momenta, 



ergodicity is straightforward to establish for the Metropolized integrator. 



Theorem 2.4 (Ergodicity). Assume 2.1 For all h > 0, observables g : O — > R, 
and initial conditions x e Q, 



/" 

where X = x. 



lim i I g(X [t/hi )dt= [ g(x)^(dx). 

_s> °° -I Jo ' ' Jn 



Proof. This proof is terse. We prove this statement for the Metropolized integra- 
tor based on a trivial partition of the molecular system. For more details please 
see [9 1 and references therein. The Metropolized integrator by construction sat- 
isfies property (PI). Moreover, its acceptance probability is strictly less than one 
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everywhere, and the integrator without accept or reject steps admits a smooth tran- 
sition density when sampled every other step. These two observations imply that 
the smooth part of the two-step transition probability of the Metropolized integra- 
tor is supported everywhere, and hence, the chain is irreducible. Irreducibility and 
property (PI) together imply ergodicity p3|. □ 



Theorems 2.3 and 2.4 are sufficient to establish that the Metropolized integra- 



tor can estimate dynamics along equilibrium trajectories of ( |T8| ). For example, as 
a corollary to the above theorems one can prove that the Metropolized integrator 
can be used to compute equilibrium correlation functions. 



Corollary 2.5. Assume \2J\ For every T > there exist positive constants h c and 
C (T) such that 

\A h (r) - A(h[r/h\)\ < C(T)h, 
for all h < h c and r G [0, T]. 

A proof of this result is provided in the Appendix. 

2.4 Case of Multiple-Time-Stepsizes 

Now we present an implementation of the patch to an explicit integrator for ([3]) 



based on a multiple-time-stepsize integrator known as RESPA [28 48 1. RESPA 
was proposed for molecular systems at constant energy in p9| . This integrator 
is a version of Verlet adapted to molecular systems with multiple time scales. It 
is designed to overcome a time-stepsize restriction imposed by rapidly changing 
short-range interactions. The scheme evaluates short-range forces on a smaller 
time-stepsize, and hence, more frequently than long-range forces. The overall 
accuracy of the algorithm is dictated by the largest time-stepsize. For Hamiltonian 



ODEs RESPA is known to exhibit resonance instabilities as described in JT7J. The 
resonance occurs between forces evaluated at the coarse time-stepsize and the 
normal modes of the molecular system excluding long-range interactions. These 
numerical instabilities also appear in generalizations of RESPA to Langevin SDEs 



[16 23 31 1. In the following we show how to patch such a generalization to tackle 
such instabilities. 

Again we partition the molecular system into m sets of particles so that each 
set has v± degrees of freedom and Y^iLi v % = v - We assume that the potential 
energy can be decomposed into fast and slow interactions: 

U(q) = U slow (q) + U fast (q) . 
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We further assume fast and slow potential forces are respectively evaluated at hf 
and h time increments. The small time-stepsize is typically a fraction of the large 



time- step size. The integrator ( [13] ) will be used, except that 9hj is replaced with 
a Metropolized version of RESPA that separately updates each element of the 
partition. Given a small time-stepsize hf, large time-stepsize h, and input 

Xo = ((Ql,--- ,Qm),(Pl,--- ,Pm)) ■ 

Set Nf — [h/hf\. The map 9hj computes a proposed move 



as 



a - {(Qi,--- ,qj-i,q*,q j+ i, 
using a RESPA update 



(P 



'Pj-liPjiPj+lJ ' ' ' iPm)) ' 



Pi.± = Pj - |Vo 3 .CW(9l, • • • , 9r 



1i 



N f 



Pj* - ifV Q .[/ fast (<Z 



1; 



g, s+i = q -■ fc + fe/M, 1 P - ft+i , 



p - 



J,JV / 



y VQ i C/ fast (q 1 , • • • , g^j, q jt k±i,q j+1 , ••■ ,q r 



1j = 9jA 

P*j = Pja ~ ^Qi U s\ 0W (qi, ■■■ , Qj-i,Qj, Qj+i, ■ ■ ■ , q m ) , 

where the inner index k runs from 1, • • • , Nf — 1 (where Nf is the number of steps 
taken at the small time-stepsize hf) and the outer index j runs from 1, • • • , m 
(where m is the number of elements in the partition). This proposed move is 
accepted with probability 



lAexpH3[tf(x*)-tf(a: )]) 



(19) 



If this proposal is rejected the momentum is reversed as before. This generaliza- 
tion of RESPA achieves a speedup when the slow potential force is expensive to 
evaluate relative to the fast potential force. 



2.5 Case of Holonomic Constraints 

Here we show how to implement the patch to molecular systems with holonomic 
constraints. Similar issues are treated in pTj[27[ from the viewpoint of sampling 
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in the presence of constraints. Again we partition the molecular system with mass 
matrix M into m sets of particles so that each set has Ui degrees of freedom and 
Y^iLi u i — v - The main difference to the previous cases is the possibility that the 
dynamics along each element of this partition is not well-defined. For instance, 
if the constraint couples all atoms in the molecule, then only the trivial partition 
will lead to well-posed dynamics. Or, if the constraint couples pairs of molecules, 
then only partitions in terms of these pairs are permissible. 

To rule out this possibility, we assume that the jth set of the partition has an 
associated scalar constraint function gi : T u — > R independent from the positions 
of the other sets. The case of vectorial constraints can be handled quite simi- 
larly, but for clarity we will consider scalar constraints here. The intersection of 
the zero level sets of these constraint functions defines the constraint manifold 
E = fl™ =1 (7^ 1 (0) C fi. The velocities of the constrained atoms are tangent to 
this manifold. These observations motivate introducing the set of all constrained 
positions and momenta denoted by TE which is known as the cotangent manifold: 

TE = {(q,p)eQ | ft(g) = 0, Vq^qfM^p, = 0, z = l,---,m}, 

where Mi is the subset of the global mass matrix M associated to the ith set of 
particles. 

For a molecular system with constraints, the probability distribution ([2]) is 
replaced by: 

dix TS (q,p)^Z- 1 e-P H ^da T x(q,p), Z = [ e^ H ^da TS (q,p) . (20) 

Here the measure cr^s represents standard volume measure on the manifold TE. 
Introduce the Lagrange multiplier Aj(i) G R. The stochastic dynamics of the 
constrained molecular system is assumed to be of the form: 

dQi =M7 ip 

dt 

dPi = -V Qi U(Q)dt + V Qi9l (Q)d^ + drtiiQ, P) > (21) 

9i(Q) =0, 



where i enumerates the elements in the partition. To check that p0| ) is an invariant 



measure of (21), its suffices to show its density is a stationary solution of the 



corresponding Fokker-Planck equation. We will assume the solution to plj ) is 
ergodic with respect to the probability distribution To eliminate the Lagrange 
multiplier appearing in pT) , 'differentiate' the constraint function twice along a 
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path of (Q,P) following either the heuristic approach in [50] or the rigorous 
approach described in [|7J 14 1. 

We will now introduce a constrained version of ( [TTj ). It is obtained by splitting 
pTj ) into a constrained Hamiltonian system: 

rdQi 



dt 
dPj 

dt 

U(Q) 



M7 1 P i 



-V Qi C/(Q)+V Q<ft (Q)A 



(22) 



0, 



and constrained thermostat dynamics: 

= 0, 




drji + V Qi gi(Q)dX it2 
V Qi9i {Q) T Mi l P t . 



(23) 



We approximate the solution to ( |22| ) by using a constrained version of Verlet 
known as RATTLE (40J. As before to maintain scalability we will separately 
propagate each set of the partition using RATTLE. Since RATTLE moves are 
time-reversible and volume-preserving (26J, a Metropolis method based on a RAT- 
TLE proposed move and the probability distribution ( f20| ) yields an acceptance 
probability that is a function of the change in energy induced by the separate RAT- 
TLE moves. As before we compose this map with an approximation to ([23]) which 
we denote by 4>t k+1 ,t k '■ TT, — > TS. The step-by-step procedure to implement this 
algorithm is given below. 

If we assume that the thermostat in ( |2"3] ) is given by Langevin dynamics, then 
its pathwise unique flow is given by: 



in 



■ {(«<,P<)}£i ^ {{Qi,e 



where we have introduced the random vector: 

"tk+l 



Pi + 4 k+u t k ))Yi 



(24) 



77 



7 P i (q)M- 1 (i feH 



i(q)dWi(s) , 



tk 



and the projection matrix: 

P t (q) = I-{V Qi gd[{V Ql 9d T MT\V Qi gdV{V Ql 9 l ) T M- 1 . 



Here / is the z/j x i/j identity matrix. To derive ( [24] ) eliminate the Lagrange multi- 
plier in ([23]) by differentiating the momentum constraint and using as an integrat- 
ing factor the matrix exp(7Pj(qr)M"^ 1 t). 
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Remark 2.6. When the mass matrix is not the identity, ( [24] ) requires computing the 
exponential of a position dependent matrix. This calculation is non-trivial, and can 
be avoided by using SHAKE in place of RATTLE and an unconstrained Ornstein- 



Uhlenbeck update in place of ( |24| ). The resulting Metropolized integrator will 
satisfy property (PI) with respect to a constrained equilibrium distribution, but 
with unconstrained velocities. It will also satisfy property (P2). 

Given a time-stepsize h, initial condition x E £1, and uniform random num- 
bers Q ~ [7(0, 1) for j = 1, • • • ,m, one step of the Metropolized integrator 



determines X\ using ( [13] ), but with the proposed moves in 9hj obtained by RAT- 
TLE as follows. Introduce the discrete Lagrange multipliers A^i, \j t 2 £ K. The 
integrator inputs a point in phase space on the constraint manifold 

x = ((<Zi,-- - ,<Z TO )>(Pi>-- - , 
and outputs a proposed move on the constraint manifold 

,Pm)) » 





= ((9!," 


• ,qj-i,qj,qj+i,- ■ 


■ ,9m), (Pi,' ■ ■ ,Pj 


-l'Pj'Pj+1' ' ' 
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(25) 



The Lagrange multiplier A^i enforces that the jth proposed positions satisfy the 
position constraint, and Xj >2 enforces that the jth proposed velocities are tangent 
to the constraint manifold. This proposed move is accepted with probability 

1 Aexp(-P[H(x*)-H(x )]) . (26) 

This acceptance probability is a function of the change in energy induced by the 
RATTLE proposed move. If this proposal is rejected the velocity is reversed, but 
remains tangent to the constraint manifold. 
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3 Applications 



This section tests the Metropolized integrators introduced in ^2} 

3.1 Lennard- Jones Fluid 

A Lennard- Jones fluid consists of n identical particles with pairwise interactions 
given by a Lennard- Jones potential energy. In what follows we use dimensionless 
units to describe this system. In these units mass is rescaled by the mass of an 
individual particle (so that the particles have unit mass), energy by the depth of 
the Lennard- Jones potential energy, and length by the point where the potential 



energy is zero. We follow the notation and setup provided in Part I of [ 1 8 1 . 

We simulate the Lennard- Jones fluid in a fixed periodic box which we call the 
simulation box. We used a truncated version of the Lennard- Jones potential in 
which the energy between two particles a distance r apart is kept constant after a 
certain cutoff distance and is given by: 

UUr) = \ S n (r) - m T t T ". 

I otherwise , 

where we have introduced /(r) = 4(l/r 12 — 1/r 6 ) and r c is bounded above by the 
size of the simulation box. Other shifts can be used to make the higher derivatives 



of Ulj continuous. The error introduced by the truncation in ( |27| ) is proportional 
to the density of the molecular system and can be made arbitrarily small by se- 
lecting the cutoff distance sufficiently large. 

The pair potentials are a function of the distance between the ith and jth parti- 
cle. If the position of the ith and jth particles are and q,, and the length of the 
simulation box i, then this distance is given by: 

ri,j(q) = KQi-Qj) mod i\ . 

In terms of this pairwise distance, the potential energy of a Lennard- Jones fluid is 
a sum of interactions between all pairs of particles: 

n— 1 n 

U ^) = H E UlMM)- (28) 

i=l j=i+l 

Evaluating the potential force requires 0(n 2 ) operations (where v is the dimension 
of configuration space), and typically dominates total computation cost. 
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3.2 Autocorrelation of Lennard- Jones Fluid 



Here we test the accuracy of the Metropolized integrator ( [13] ) based on trivial and 
per particle partitions of the molecular system. In the former proposed moves in 
the Metropolis steps are obtained by per particle Verlet updates, and in the latter 
by global Verlet updates. The numerics indicate that the order of accuracy of both 
methods is roughly 0(h 2 ) and that the error constant of the per particle partition 
is approximately an order of magnitude smaller. 

To test accuracy we use the Metropolized integrator to estimate the equilib- 
rium momentum autocorrelation of the Lennard- Jones fluid. For the numerical 
experiment, we set the fluid's density g = 0.8442 and temperature T = 0.728 
(units are dimensionless as described earlier). For these values the phase of the 
Lennard- Jones fluid is liquid, and close to the triple (gas-liquid- solid) point. We 
also fix the number of particles to be n = 25 and the degrees of freedom to be 
d = 2. The size of the simulation box is i = (n/ q) 1 / 2 « 5.04. A reasonable cutoff 
distance in ( [27] ) at the selected density is r c = 2.5. The initial positions of the par- 



ticles are chosen to be the vertices of a square lattice that fills the simulation box. 
For instance, the length of each square in the lattice can be chosen to be 1/ [n 1 / 2 ] . 
The initial velocities are sampled from the Maxwell distribution. The thermostat 
parameter is set equal to 7 = 1.0. 

In the simulations we estimate the true velocity autocorrelation over a time- 
interval [0, 1]. Set N = [T/h\ and introduce an evenly-spaced mesh in computa- 
tional time t k = hk for all < k < N. Let A h : [0, 1] -> R denote the velocity 



correlation function obtained by the Metropolized integrator ( 13 1: 



1 N 

A " (r) = iToo N £ P W/>j P L(**)/M > (r > 0) . 



fc=i 



Define the relative Richardson error as 

def SU Pre[0,l] 



A h (r) - A 2h (r) 



(29) 



su Pre[o,i] l A ( r )l 
An empirical estimate of is plotted for time-stepsizes 

h = {0.005,0.0025,0.00125,0.000625} 

in Fig.[T]with N = 10 8 . The denominator in e h is calculated by using the approxi- 
mation of A at h = 0.000625. The figure shows that the error of the Metropolized 
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integrator is approximately 0(h 2 ). Moreover it shows that the Metropolized in- 
tegrator based on a per particle partition is nearly an order of magnitude more 
accurate than the algorithm based on a trivial partition. 



3.3 Scaling of Metropolized Integrator 

Until now the numerics dealt with a fixed number of particles in d = 2 dimensions. 



Next we address how the integrator ( [13] ) scales with the number of particles keep- 
ing the density and temperature fixed. As before we consider two types of parti- 
tions: per particle and trivial. This time we will assume the particles are in three- 
dimensional space. Recall that in the former proposed moves in the Metropolis 
steps are obtained by per particle Verlet updates, and in the latter by global Verlet 
updates. We will show that the type of proposed move affects the scalability of 
the Metropolized integrator. In particular, the per particle partition will lead to a 
scalable algorithm. 

By fixing the density and temperature, the stiffness of the molecular system 
is fixed. Thus, from the viewpoint of numerical analysis, the time-stepsize ought 
to be independent of system size. However, the acceptance probability in the 
Metropolized integrator depends on the change in energy induced by the pro- 



posal move (see (|T5])). If this proposal move is global, then the magnitude of this 
change in energy increases with system size. Hence, the acceptance probability is 
inversely related to system size, in general. This poor scaling of Metropolis meth- 



ods based on global moves is well-known in the literature (see, e.g., p} |32l[38| ). 

For the numerical experiment, we set the fluid's density g = 0.8442 and tem- 
perature T = 0.728 (as before units are dimensionless). The mean acceptance 
probability is computed along a long time trajectory of the integrators (10 6 steps) 
with a fixed time-stepsize of h = 0.01 and a variety of system sizes. The initial 
positions of the particles are the vertices of a cubic lattice contained in the sim- 
ulation box. The length of each cube is given by [(1/f?) 1 / 3 ] (g = 0.8442 is the 
density of the fluid). The initial velocities are sampled from the Maxwell distribu- 
tion. The thermostat parameter and Lennard- Jones cutoff distances are set equal 
to 7 = 1.0 and r c = 2.5, respectively. 

Figure [2] shows the outcome of the experiment: the mean acceptance prob- 
ability per particle for the per particle and trivial partitions as a function of the 
number of particles. The acceptance probability for the trivial partition clearly 
deteriorates with system size. On the other hand, the acceptance probability for 
the per particle partition is independent of system size. In fact, the acceptance 
probability per particle is equal to one to within round off error. This result seems 
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to defy intuition since the particles experience Lennard- Jones interaction. But, 
keep in mind that these interactions are, in fact, short-range due to the cutoff in 
the Lennard- Jones potential energy. Hence, the change in energy induced by per 
particle Verlet moves is independent of system size. 



3.4 Autocorrelation of Lennard- Jones Dumbbells 

A rigid dumbbell is a type of molecule that involves a holonomic constraint. It 
consists of a pair of particles constrained to a fixed distance from one another. The 
constraint arises when the spring joining a flexible dumbbell infinitely stiffens. 
Consider n identical dumbbells where particles in separate dumbbells interact via 
a Lennard- Jones potential energy ( |28| ). 

For the numerical experiment, we simulate the dumbbells in a simulation box 
of length I and d = 2 dimensions as before. We set the system's density g = 0.998 
and temperature T = 3.0 (units are dimensionless as described earlier). We also 
fix the number of dumbbells to n = 30 and the length of each dumbbell to £ = 1. 
If the positions of the ith pair of particles describing the ith dumbbell are q i 1 and 
q i 2 , then the constraint function associated to the ith dumbbell is given by: 

9i(q) = \(Qi,i -Qi,2) mod£| 2 -£o, i = 1, • • • ,n . 
The size of the simulation box is £ = (2n/g) 1 ^ 2 ps 7. The cutoff distance in ( |27| ) 



is set at r c = 3.0. The initial positions of the dumbbells are chosen randomly in 
the simulation box, but with no overlap. The initial velocities are sampled from 
the Maxwell distribution constrained to be tangent to the constraint manifold at 
the initial positions and temperature. The thermostat parameter is set to 7 = 1.0. 

To quantify the accuracy of the numerical method, we use the Metropolized 
integrator based on a per dumbbell partition. We again use the relative Richardson 
error defined in ( [29] ) to estimate the rate of convergence of the method. Figure [3] 



shows the velocity autocorrelation that we wish to approximate over the time- 
interval [0, 1]. It is noticeably different from the velocity autocorrelation in the 
case of the Lennard- Jones cluster without constraints. The figure also shows a plot 
of the Richardson error as a function of time-stepsize. The rate of convergence 
appears to be approximately 0(h 2 ). 
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4 Conclusions 



This paper showed how thinking probabilistically helps to design good integrators 
for SDEs arising in molecular dynamics. The paper brought ideas from Monte 
Carlo into molecular dynamics, to obtain a Metropolized integrator that is uncon- 
ditionally stable, pathwise accurate, and still explicit. While the examples and 
theory in the paper focused on Langevin dynamics, the methodology can be ex- 
tended to other thermostats. The patch we propose is simple and, as we showed, 
its computational overhead is minimum compared to an 'unpatched' integrator. 

An open question about the Metropolized integrator is its convergence rate to 
equilibrium. In general, it is not expected that the Metropolized integrator inherits 
all of the mixing properties of the exact solution to the SDE. The main reason 
is conditional stability of the underlying Verlet integrator that generates proposal 
moves. Indeed for any time-stepsize one can find an energy value above which the 
drift in this integrator gives proposed moves that increase the energy, in contrast 
to the exact drift in the SDE which always centers the solution towards lower 
energy values. Since higher energy values have a lower equilibrium probability 
weight, these proposed moves are typically rejected. While these rejections ensure 
that the Metropolized integrator is ergodic, at high energy values they prevent the 
integrator from inheriting all of the mixing properties of the true solution. 

This issue has been addressed for the MALA algorithm in the overdamped 
limit of Langevin dynamics [6|. It turns out that MALA converges to its equilib- 
rium distribution at an exponential rate up to terms exponentially small in time- 
stepsize. However, in the overdamped limit positions are no longer differentiable, 
momentum is no longer present, and the MALA algorithm does not involve mo- 
mentum flips. Intuitively one expects these momentum flips to reduce stagnation 
at high energy values, and hence, enhance the mixing rate. Future research will 
investigate the effect of these momentum flips to this mixing rate. 

Another open question concerns the acceptance probability of the Metropolized 
integrator as the dimension tends to oo. The numerical experiment in Figure [2] in- 
dicates that the acceptance probability of the Metropolized integrator based on 
global Verlet moves scales as (9(n -1 / 6 ) where n is the number of particles. Thus, 
to obtain an 0(1) acceptance probability the time stepsize should be proportional 
to n -1 / 6 , and the integrator would require O^n 1 ^ 6 ) steps to traverse phase space. 
This scaling property would imply that the Metropolized integrator is more ef- 
ficient at making 0(1) moves in state space as compared with the (9(n 1//4 ) and 
0(n 1 ^ 3 ) scalings of respectively hybrid Monte-Carlo and MALA [|5j. This ques- 
tion will also be investigated in future work. 
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Appendix 

Proof of Corollary 2.5 We prove this estimate for the equilibrium correlation of 
any Lipschitz function g £ C 0,1 (f2,M). This assumption includes the special 
case of a scalar velocity autocorrelation function. Ergodicity of the Metropolized 
integrator (see Theorem [24] ) implies that 

1 N 

j^^S^ x Lv^jM x L(* fc +r)//>j) = (g(Xit/h\)g(x l(t+T)/hl )) 

k=l 

for any t > 0. (Recall, that the angle brackets denote a double average with 
respect to initial conditions distributed according to the equlibrium distribution of 
< p~8] ) and realizations of the Wiener process.) Thus, we wish to estimate: 

e h ^ \{g(X mi )g(X l{t+T)/hi )) - (g(Y(t))g(Y(t + h[r/h\)))\ 
= \(g(x)(g(X lT/hl )-g(Y(h[r/h\))))\. 

By the Cauchy-Schwarz inequality, 

e h < C a (j \E x {g(X [T/hl ) -g(Y(h[r/h\))}\ V ' 

where C g = (J n \g(x)\ 2 dfi) 1 ^ 2 . By Jensen's inequality, 

e h < C g (w{\g(X lT/hi ) -g(Y(h[r/h\))\ 2 } 
The Lipschitz assumption on g implies, 

\g(x) -g(y)\ < L g \x-y\ 



1/2 
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for some constant L g > and for every x, y G Vt. Hence, 

1/2 

e h < L g C g (E»{\X [T/hi - Y(h[r/h\)\ 2 y 



Pathwise convergence of the Metropolized integrator (see Theorem 2.3 ) implies 

1 /2 

(E" { | X lT/h} -Y(h[r/h\)\ 2 ^ <C{T)h. 

From which it follows, that the accuracy of the Metropolized integrator in com- 
puting the equilibrium correlation in g is 0(h) with a prefactor that depends on 
the function g. □ 

References 

[1] E. Akhmatskaya and S. Reich, GSHMC: An efficient method for molecular 
simulation, J. Comput. Phys. 227 (2008), 4937-4954. 

[2] M. P. Allen and D. J. Tildesley, Computer simulation of liquids, Clarendon 
Press, 1987. 

[3] H. C. Andersen, Molecular dynamics simulations at constant pressure 
and/or temperature, J. Chem. Phys. 72 (1980), 2384. 

[4] G. Benetin and A. Giorgilli, On the Hamiltonian interpolation of near to 
the identity symplectic mappings with applications to symplectic integration 
algorithms, J. Statist. Phys. 74 (1994), 1117-1143. 

[5] A. Beskos, N. S. Pillai, G. O. Roberts, J. M. Sanz-Serna, and A. M. Stuart, 
Optimal tuning of hybrid Monte-Carlo, arXiv: 1001.4460, 2010. 

[6] N. Bou-Rabee, M. Hairer, and E. Vanden-Eijnden, Non-asymptotic mixing 
of the MALA algorithm, arXiv: 1008.35 14vl [math.PR], 2010. 

[7] N. Bou-Rabee and H. Owhadi, Stochastic variational integrators, IMA J. of 
Numer. Anal. 29 (2009), 421-443. 

[8] N. Bou-Rabee and E. Vanden-Eijnden, Reconciling Monte Carlo and molec- 
ular dynamics, Preprint, 2009. 



22 



[9] N. Bou-Rabee and E. Vanden-Eijnden, Pathwise accuracy and ergodicity of 
Metropolized integrators for SDEs, CPAM 63 (2010), 655-696. 

[10] A. Briinger, C. L. Brooks, and M. Karplus, Stochastic boundary conditions 
for molecular dynamics simulations of ST2 water, Chem. Phys. Lett. 105 
(1984), 495-500. 

[11] G. Bussi, D. Donadio, and M. Parrinello, Canonical sampling through ve- 
locity rescaling, J. Chem. Phys. 126 (2007), 014101. 

[12] G. Bussi and M. Parrinello, Stochastic thermostats: Comparison of local and 
global schemes, Computer Physics Communications 179 (2008), 26-29. 

[13] E. Cances, F. Legoll, and G. Stoltz, Theoretical and numerical comparison 
of some sampling methods for molecular dynamics, Mathematical Modelling 
and Numerical Analysis 41 (2007), 351-389. 

[14] G. Ciccotti, T. Lelievre, and E. Vanden-Eijnden, Projections of diffusions on 
submanifolds: Application to mean force computation, CPAM 61 (2008), 
0001-0039. 

[15] W. E and D. Li, The Andersen thermostat in molecular dynamics, CPAM 61 
(2008), 96-136. 

[16] W. Fong, Multi-scale methods in time and space for particle simulations, 
Ph.D. thesis, Stanford University, 2009. 

[17] W. Fong, E. Darve, and A. Lew, Stability of asynchronous variational inte- 
grators, J. Comput. Phys. 227 (2008), 8367-8394. 

[18] D. Frenkel and B. Smit, Understanding molecular simulation: From algo- 
rithms to applications, second edition, Academic Press, 2002. 

[19] J. Goodman and A. Sokal, Multigrid Monte Carlo method: Conceptual foun- 
dations, Phys. Rev. D 40 (1989), 2035-2071. 

[20] E. Hairer, C. Lubich, and G. Wanner, Geometric numerical integration, 
Springer Series in Computational Mathematics, vol. 31, Springer, 2006. 

[21] C. Hartmann, An ergodic sampling scheme for constrained Hamiltonian sys- 
tems with applications to molecular dynamics, J. Stat. Phys. 130 (2008), 
687-712. 



23 



[22] D. J. Higham, X. Mao, and A. M. Stuart, Strong convergence of Euler-type 
methods for nonlinear stochastic differential equations, SIAM J. Numer. 
Anal. 40 (2002), 1041-1063. 

[23] J. A. Izaguirre, D. P. Catarello, J. M. Wozniak, and R. D. Skeel, Langevin 
stabilization of molecular dynamics, J. Chem. Phys. 114 (2001), 2090. 

[24] B. Leimkuhler, E. Noorizadeh, and F. Theil, A gentle stochastic thermostat 
for molecular dynamics, J. Stat. Phys. 135 (2009), 261-277. 

[25] B. Leimkuhler and S. Reich, Simulating Hamiltonian dynamics, Cambridge 
Monographs on Applied and Computational Mathematics, Cambridge Uni- 
versity Press, 2004. 

[26] B. Leimkuhler and R. Skeel, Symplectic numerical integrators in con- 
strained Hamiltonian systems, JCP 112 (1994), 1 17-125. 

[27] T. Lelievre, M. Rousset, and G. Stoltz, Langevin dynamics with constraints 
and computation of free energy differencs, arXiv:1006.4914vl, 2010. 

[28] A. Lew, J. E. Marsden, M. Ortiz, and M. West, Asynchronous variational 
integrators, Arch. Ration. Mech. An. 167 (2003), 85-145. 

[29] D. Li, On the rate of convergence to equilibrium of the Andersen thermostat 
in molecular dynamics , J. Stat. Phys. 129 (2007), 265-287. 

[30] J. S. Liu and Y. N. Wu, Generalised Gibbs sampler and multigrid Monte 
Carlo for Bayesian computation, Biometrika 87 (2000), 353-369. 

[31] Q. Ma, J. A. Izaguirre, and R. D. Skeel, VERLET-I/R-RESPA/IMPULSE is 
limited by nonlinear instabilities, SIAM J. Sci. Comput. 24 (2003), 1951- 
1973. 

[32] J. C. Mattingly, N. S. Pillai, and A. M. Stuart, Diffusion limits of the random 
walk Metropolis algorithm in high dimensions, arXiv: 1003.4306, 2010. 

[33] K. L. Mengersen and R. L. Tweedie, Rates of convergence of the Hastings 
and Metropolis algorithms, Ann. Stat. 24 (1996), 101-121. 

[34] N. Metropolis, A. W. Rosenbluth, A. H. Teller, and E. Teller, Equations of 
state calculations by fast computing machines, J. Chem. Phys. 21 (1953), 
1087-1092. 



24 



[35] G. N. Milstein and M. V. Tretyakov, Numerical integration of stochastic dif- 
ferential equations with nonglobally Lipschitz coefficients, SIAM J. Numer. 
Anal. 43(2005), 1139-1154. 

[36] J. Moser, Lectures on Hamiltonian systems, vol. 81, Mem. AMS, 1968. 

[37] S. Reich, Backward error analysis for numerical integrators, SIAM J. Num. 
Anal. 36 (1999), 1549-1570. ' 

[38] G. O. Roberts and J. S. Rosenthal, Optimal scaling of discrete approxima- 
tions to Langevin diffusions, J. Roy. Statist. Soc. Ser. B 60 (1998), 255-268. 

[39] G. O. Roberts and R. L. Tweedie, Geometric convergence and central 
limit theorems for multidimensional Hastings and Metropolis algorithms, 
Biometrika 1 (1996), 95-110. 

[40] J. Ryckaert, G. Ciccotti, and H. Berendsen, Numerical integration of the 
Cartesian equations of motion of a system with constraints: Molecular dy- 
namics of n-alkanes, JCP 23 (1977), 327-341. 

[41] A. A. Samoletov, M. A. Chaplain, and C. P. Dettmann, Thermostats for 
"slow" configurational modes, J. Stat. Phys. 128 (2008), 1321-1336. 

[42] A. Scemama, T. Lelievre, G. Stoltz, E. Cances, and M. Caffarel, An efficient 
sampling algorithm for variational Monte Carlo, J. Chem. Phys. 125 (2006), 
114105. 

[43] T. Schneider and E. Stoll, Molecular-dynamics study of a three-dimensional 
one-component model for distortive phase transitions, Physical Review B 17 
(1978), 1302-1322. 

[44] G. Stoltz, Some mathematical methods for molecular and multiscale simu- 
lation, Ph.D. thesis, Ecole Nationale des Ponts et Chaussees, 2007. 

[45] D. Talay, Simulation and numerical analysis of stochastic differential sys- 
tems : a review, Probabilistic Methods in Applied Physics (P. Kree and 
W. Wedig, eds.), vol. 451, Springer- Verlag, Berlin, 1995, pp. 54-96. 

[46] , Stochastic Hamiltonian systems: Exponential convergence to the 

invariant measure, and discretization by the implicit Euler scheme, Markov 
Processes and Related Fields 8 (2002), 1-36. 



25 



[47] D. Talay and L. Tubaro, Expansion of the global error for numerical schemes 
solving stochastic differential equations, Stoch. Anal. Appl. 8 (1990), 94- 
120. 

[48] M. E. Tuckerman and B. J. Berne, Stochastic molecular dynamics in systems 
with multiple time scales and memory friction, J. Chem. Phys. 95 (1991), 
4389-4396. 

[49] M. E. Tuckerman, B. J. Berne, and G. Martyna, Reversible multiple time 
scale molecular dynamics, J. Chem. Phys. 97 (1992), 1990-2001. 

[50] E. Vanden-Eijnden and G. Ciccotti, Second-order integrators for Langevin 
equations with holonomic constraints, Chem. Phys. Lett. 429 (2006), 310- 
316. 

[51] E. Vanden-Eijnden and M. Venturoli, Markovian milestoning with Voronoi 
tessellations, J. Chem. Phys. 130 (2008), 194101. 

[52] , Exact rate calculations by trajectory parallelization and twisting, 

In press, 2009. 

[53] L. Verlet, Computer "experiments" on classical fluids. I. Thermodynamical 
properties of Lennard- Jones molecules, Phys. Rev. 159 (1967), 98-103. 



26 




o 

LU 

c 

o 

CO 
T3 

CO 
sz 
g 

DC 

CD 
> 

CD 

DC 



10" 



10" 



10" 









— Trivial 




— ©— Per-Particle 




--h 2 



10" 



Figure 1 : Autocorrelation of Lennard- Jones Fluid. The panel on the top shows the 
velocity autocorrelation function over the time interval [0, 1]. The panel on the bottom 



shows a loglog plot of the Richardson error ((|29[)) of the Metropolized integrator based on 
trivial and per-particle partitions. The molecular system integrated is the Lennard- Jones 
fluid with 25 particles in a square box. The dashed line represents a reference slope of 
h 2 . A total of N = 10 s samples were used to obtain these empirical estimates. Notice 
that the Metropolized integrator based on per-particle Verlet moves is about an order of 
magnitude more accurate than the Metropolized integrator based on global Verlet moves. 
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Figure 2: Mean Acceptance Probability as a Function of System Size. This figure 
shows the mean acceptance probability per particle for the integrator ( fT3] l based on a 
trivial and per particle partition. From the plot it is clear that a) the acceptance probability 
for the integrator based on a per particle partition is independent of system size and b) the 
acceptance probability for the trivial partition scales like n^ 1 / 6 where n is the number of 
particles. 
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Figure 3: Autocorrelation of Lennard- Jones Dumbbells. The panel on the top shows 
the velocity autocorrelation function. The panel on the bottom shows a loglog plot of an 
empirical estimate of the Richardson error (((29])) of the Metropolized integrator based on 
per-dumbbell RATTLE proposed moves. The molecular system integrated is the Lennard- 
Jones system of 15 dumbbells in a square box. The dashed line represents a reference 
slope of h 2 . A total of N = 10 8 samples were used to obtain this estimate. 
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